library_load <- suppressMessages(
list(
# Seurat
library(Seurat),
# Scran doubletFinder implementation
library(scDblFinder),
# FastMNN implementation
library(SeuratWrappers),
# Data
library(tidyverse),
# Reticulate
library(reticulate),
# Biomart
library(biomaRt),
# Plotting
library(ggplot2),
library(knitr),
library(kableExtra),
library(IRdisplay)
)
)
random_seed <- 42
set.seed(random_seed)
options(warn=-1)
# Set working directory to project root
setwd("/research/peer/fdeckert/FD20200109SPLENO")
# Source files
source("plotting_global.R")
source("bin/seurat_qc.R")
source("bin/seurat_dea.R")
# Files
so_raw_file <- "data/object/raw.rds"
so_qc_rds <- "data/object/qc.rds"
so_qc_h5ad <- "data/object/qc.h5ad"
# Plotting Theme
ggplot2::theme_set(theme_global_set()) # From project global source()
so_raw <- readRDS(so_raw_file)
so_raw$tissue <- factor(so_raw$tissue, levels=c("Myeloid", "Progenitor"))
so_raw$treatment <- factor(so_raw$treatment, levels=c("NaCl", "CpG"))
GetAssayData(so_raw, assay="RNA", slot="data")["Bmp4", ] %>% sum
Empty droplets were determined with CellRanger V3.0.2 Lun et al., 2019 EmptyDrop heuristic. RNAse activity of granulocytes might be wrongly identified as empty cells by CellRanger.
Typical Sample A steep drop-off is indicative of good separation between the cell-associated barcodes and the barcodes associated with empty GEMs. A ideal barcode rank plot has a distincitve shape, which is referred to as a "cliff and knee".
Heterogeneous Sample Heterogeneous populations of cells in a sample result in two "cliff and knee" distributions. However, there should still be clear separation between the bacodes.
Compromised Sample Round curve and lack of steep cliff may indicate low sample quality or loss of single-cell behavior. This can be due to a wetting failure, premature cell lysis, or low cell viability.
Compromised Sample Defined cliff and knee, but the total number of barcodes detected may be lower than expected. This can be caused by a sample clog or inaccurate cell count.
options(repr.plot.width=10, repr.plot.height=5)
rank_plot_qc_1 <- rank_plot_qc(so_raw, color_color=color$cellranger_class, formular=tissue~sample_rep+treatment)
rank_plot_qc_1
so_qc <- subset(so_raw, subset=cellranger_class=="Cell")
# Get Seurat cell cycle genes
cc_genes_seurat_s <- cc.genes.updated.2019$s.genes
cc_genes_seurat_g2m <- cc.genes.updated.2019$g2m.genes
# Get mouse orthologs from human gene simbols
httr::set_config(httr::config(ssl_verifypeer=FALSE))
hgnc_mart <- useMart("ensembl", dataset="hsapiens_gene_ensembl", host="https://dec2021.archive.ensembl.org/")
mm_mart <- useMart("ensembl", dataset="mmusculus_gene_ensembl", host="https://dec2021.archive.ensembl.org/")
cc_genes_seurat_s <- getLDS(attributes=c("hgnc_symbol"), filters="hgnc_symbol", values=cc_genes_seurat_s, mart=hgnc_mart, attributesL=c("mgi_symbol"), martL=mm_mart, uniqueRows=TRUE)
cc_genes_seurat_s <- cc_genes_seurat_s[, 2]
cc_genes_seurat_g2m <- getLDS(attributes=c("hgnc_symbol"), filters="hgnc_symbol", values=cc_genes_seurat_g2m, mart=hgnc_mart, attributesL=c("mgi_symbol"), martL=mm_mart, uniqueRows=TRUE)
cc_genes_seurat_g2m <- cc_genes_seurat_g2m[, 2]
# Save cell cycle genes
saveRDS(list(s=cc_genes_seurat_s, g2m=cc_genes_seurat_g2m), "data/annotation/cell_cycle_seurat/cc_genes.rds")
# Compute cell cycle score
so_qc <- CellCycleScoring(so_qc, s.features=cc_genes_seurat_s, g2m.features=cc_genes_seurat_g2m, set.ident=FALSE)
colnames(so_qc@meta.data) <- gsub("Phase", "cc_phase_class", colnames(so_qc@meta.data))
colnames(so_qc@meta.data) <- gsub("S.Score", "msS_RNA", colnames(so_qc@meta.data))
colnames(so_qc@meta.data) <- gsub("G2M.Score", "msG2M_RNA", colnames(so_qc@meta.data))
so_qc$msCC_diff_RNA <- so_qc$msS_RNA - so_qc$msG2M_RNA
so_qc$cc_phase_class <- factor(so_qc$cc_phase_class, levels=names(color$cc_phase_class))
SingleR identifies marker genes from the reference and uses them to compute assignment score (based on the Spearman correlation across markers) for each cell in the test dataset against each label in the reference. The label with the highest doublet_stat is the assigned to the test cell, possibly with further fine-tuning to resolve closely related labels.
first.labels: Labels before fine-tuning
labels: Labels after fine-tuning
pruning: labels after pruning
if(FALSE) {
# Load reference data
ref <- ImmGenData()
# Seurat object to SingleCellExperiment
sce <- SingleCellExperiment(list(counts=so_qc@assays$RNA@counts))
# # Predict labels
label_main <- SingleR::SingleR(test=sce, ref=ref, labels=ref$label.main, assay.type.test="counts", de.method="classic")
label_fine <- SingleR::SingleR(test=sce, ref=ref, labels=ref$label.fine, assay.type.test="counts", de.method="classic")
saveRDS(label_main, "result/singler/label_main_immgen.rds")
saveRDS(label_fine, "result/singler/label_fine_immgen.rds")
} else {
label_main <- readRDS("result/singler/label_main_immgen.rds")
label_fine <- readRDS("result/singler/label_fine_immgen.rds")
}
# Add labels to Seurat object
label_main_meta <- as.data.frame(label_main) %>% dplyr::select(pruned.labels, tuning.scores.first) %>% dplyr::rename(label_main_immgen=pruned.labels, delta_score_main_immgen=tuning.scores.first)
label_fine_meta <- as.data.frame(label_fine) %>% dplyr::select(pruned.labels, tuning.scores.first) %>% dplyr::rename(label_fine_immgen=pruned.labels, delta_score_fine_immgen=tuning.scores.first)
so_qc <- AddMetaData(so_qc, cbind(label_main_meta, label_fine_meta))
if(FALSE) {
# Load reference data
ref <- readRDS("data/haemosphere/se_haemosphere.rds")
# Seurat object to SingleCellExperiment
sce <- SingleCellExperiment(list(counts=so_qc@assays$RNA@counts))
# # Predict labels
label_main <- SingleR::SingleR(test=sce, ref=ref, labels=ref$label.main, assay.type.test="counts", de.method="classic")
label_fine <- SingleR::SingleR(test=sce, ref=ref, labels=ref$label.fine, assay.type.test="counts", de.method="classic")
saveRDS(label_main, "result/singler/label_main_haemosphere.rds")
saveRDS(label_fine, "result/singler/label_fine_haemosphere.rds")
} else {
label_main <- readRDS("result/singler/label_main_haemosphere.rds")
label_fine <- readRDS("result/singler/label_fine_haemosphere.rds")
}
# Add labels to Seurat object
label_main_meta <- as.data.frame(label_main) %>% dplyr::select(pruned.labels, tuning.scores.first) %>% dplyr::rename(label_main_haemosphere=pruned.labels, delta_score_main_haemosphere=tuning.scores.first)
label_fine_meta <- as.data.frame(label_fine) %>% dplyr::select(pruned.labels, tuning.scores.first) %>% dplyr::rename(label_fine_haemosphere=pruned.labels, delta_score_fine_haemosphere=tuning.scores.first)
so_qc <- AddMetaData(so_qc, cbind(label_main_meta, label_fine_meta))
I first tried the strategy from Sala et al., 2019 which was also used in Barile et al., 2021 with the source documented here. We only do the refined method with combined samples. The functions used in the script are documented here. We do doublet detection on sample groups (Myeloid and Progenitor cells of a replicated) and integrated those groups per batch and then by treatment with FastMNN. However, I default later back to just use the standard approach.
qc_class_set <- function(so) {
print(so)
# QC filter for all cells
so$pMt_RNA_min <- -Inf
so$pMt_RNA_max <- 5
so$nCount_RNA_min <- 1500
so$nCount_RNA_max <- max(so$nCount_RNA)
so$qc_class <- ifelse(
so$cellranger_class=="Cell" &
so$nCount_RNA <= so$nCount_RNA_max &
so$nCount_RNA > so$nCount_RNA_min &
so$pMt_RNA <= so$pMt_RNA_max &
so$pMt_RNA > so$pMt_RNA_min,
"pass", "fail"
)
so <- subset(so, subset=qc_class=="pass")
print(so)
return(so)
}
so_db <- Seurat::SplitObject(so_qc, split.by="sample_group")
so_db <- lapply(so_db, qc_class_set)
An object of class Seurat 14772 features across 8440 samples within 1 assay Active assay: RNA (14772 features, 0 variable features) An object of class Seurat 14772 features across 5381 samples within 1 assay Active assay: RNA (14772 features, 0 variable features) An object of class Seurat 14772 features across 11470 samples within 1 assay Active assay: RNA (14772 features, 0 variable features) An object of class Seurat 14772 features across 6730 samples within 1 assay Active assay: RNA (14772 features, 0 variable features) An object of class Seurat 14772 features across 6727 samples within 1 assay Active assay: RNA (14772 features, 0 variable features) An object of class Seurat 14772 features across 3290 samples within 1 assay Active assay: RNA (14772 features, 0 variable features) An object of class Seurat 14772 features across 10947 samples within 1 assay Active assay: RNA (14772 features, 0 variable features) An object of class Seurat 14772 features across 4567 samples within 1 assay Active assay: RNA (14772 features, 0 variable features)
# Compute doublet density per sample group
so_db <- lapply(so_db, function(so) {
DefaultAssay(so) <- "RNA"
# Remove empty genes from split
cnt <- GetAssayData(so, assay="RNA", slot="counts")
cnt <- cnt[rowSums(cnt) > 0, ]
# Compute doublet density
sce <- scDblFinder::scDblFinder(cnt)
so$doublet_score <- sce$scDblFinder.score
so$doublet_score_log2 <- log2(sce$scDblFinder.score)
so$doublet_class <- sce$scDblFinder.class
return(so)
}
)
Assuming the input to be a matrix of counts or expected counts. Creating ~5000 artificial doublets... Dimensional reduction Evaluating kNN... Training model... iter=0, 405 cells excluded from training. iter=1, 468 cells excluded from training. iter=2, 481 cells excluded from training. Threshold found:0.481 406 (7.5%) doublets called Assuming the input to be a matrix of counts or expected counts. Creating ~5000 artificial doublets... Dimensional reduction Evaluating kNN... Training model... iter=0, 604 cells excluded from training. iter=1, 493 cells excluded from training. iter=2, 448 cells excluded from training. Threshold found:0.506 396 (5.9%) doublets called Assuming the input to be a matrix of counts or expected counts. Creating ~5000 artificial doublets... Dimensional reduction Evaluating kNN... Training model... iter=0, 190 cells excluded from training. iter=1, 233 cells excluded from training. iter=2, 236 cells excluded from training. Threshold found:0.74 159 (4.8%) doublets called Assuming the input to be a matrix of counts or expected counts. Creating ~5000 artificial doublets... Dimensional reduction Evaluating kNN... Training model... iter=0, 322 cells excluded from training. iter=1, 287 cells excluded from training. iter=2, 274 cells excluded from training. Threshold found:0.504 228 (5%) doublets called
options(repr.plot.width=30, repr.plot.height=5)
# Compute sub-cluster per sample group
so_db <- lapply(so_db, function(so) {
# Cluster function
so <- ScaleData(so, features=rownames(so), assay="RNA", verbose=FALSE)
so <- FindVariableFeatures(so)
so <- RunPCA(so, npcs=50, features=VariableFeatures(so), assay="RNA", verbose=FALSE)
so <- FindNeighbors(so, dims=1:50, k.param=15, reduction="pca", verbose=FALSE)
so <- FindClusters(so, verbose=FALSE, resolution=1, algorithm=1)
so <- RunUMAP(so, dims=1:50, n.neighbors=30, verbose=FALSE)
dplot_1 <- dplot(so, reduction="umap", group_by="treatment") + scale_color_manual(values=color$treatment) + ggtitle("Treatment")
dplot_2 <- dplot(so, reduction="umap", group_by="tissue") + scale_color_manual(values=color$tissue) + ggtitle("Tissue")
dplot_3 <- dplot(so, reduction="umap", group_by="doublet_class") + ggtitle("doublet_class")
fplot_1 <- fplot(so, reduction="umap", features="doublet_score") + scale_colour_viridis_c(option="inferno") + ggtitle("doublet_score")
fplot_2 <- fplot(so, reduction="umap", features="doublet_score_log2") + scale_colour_viridis_c(option="inferno") + ggtitle("doublet_score_log2")
dplot_comb <- dplot_1 + dplot_2 + dplot_3 + fplot_1 + fplot_2 + plot_layout(ncol=5)
plot(dplot_comb)
return(so)
}
)
Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale. Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale. Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale. Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale. Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale. Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.
Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale. Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.
so_db <- merge(so_db[[1]], so_db[2:length(so_db)])
so_qc <- AddMetaData(so_qc, dplyr::select(so_db@meta.data, doublet_class, doublet_score, doublet_score_log2))
qc_class_set <- function(so) {
# QC filter for all cells
so$pMt_RNA_min <- -Inf
so$pMt_RNA_max <- 5
so$nCount_RNA_min <- 1500
so$nCount_RNA_max <- max(so$nCount_RNA)
so$qc_class <- ifelse(
so$cellranger_class == "Cell" &
so$nCount_RNA <= so$nCount_RNA_max &
so$nCount_RNA > so$nCount_RNA_min &
so$pMt_RNA <= so$pMt_RNA_max &
so$pMt_RNA > so$pMt_RNA_min &
so$doublet_class == "singlet",
"pass", "fail"
)
# Remove Ery from Myeloid sort before computing nFeature_RNA
if(so$tissue[1]=="Myeloid") {
so$qc_class <- ifelse(
so$label_main_haemosphere=="Ery" | so$label_fine_haemosphere %in% c("preCFUE", "CFUE", "pbEry", "poEry", "Retic") | so$pHb_RNA >= 5,
"fail", so$qc_class
)
}
so_tmp <- subset(so, subset=qc_class=="pass")
so$nFeature_RNA_max <- max(so_tmp$nFeature_RNA)
so$nFeature_RNA_min <- quantile(so_tmp$nFeature_RNA, 0.01)
so$qc_class <- ifelse(
so$qc_class == "pass" &
so$nFeature_RNA <= so$nFeature_RNA_max &
so$nFeature_RNA > so$nFeature_RNA_min,
"pass", "fail"
)
return(so)
}
so_qc <- Seurat::SplitObject(so_qc, split.by="sample_name")
so_qc <- lapply(so_qc, qc_class_set)
so_qc <- merge(so_qc[[1]], so_qc[2:length(so_qc)])
saveRDS(so_qc, "data/object/raw_filter_label.rds")
density_plot_qc_1 <- density_plot_qc(so=so_qc, title="Density plot UMI count", x=nCount_RNA, xlab="log10(UMI count)", min=nCount_RNA_min, max=nCount_RNA_max, fill_color=color$tissue)
density_plot_qc_2 <- density_plot_qc(so=so_qc, title="Density plot Feature count", x=nFeature_RNA, xlab="log10(Feature count)", min=nFeature_RNA_min, max=nFeature_RNA_max, fill_color=color$tissue)
density_plot_qc_3 <- density_plot_qc(so=so_qc, title="Density plot Mt %", x=pMt_RNA, xlab="Mt [%]", min=0, max=pMt_RNA_max, log10=FALSE, fill_color=color$tissue)
options(repr.plot.width=22.5, repr.plot.height=5)
density_plot_qc_1 + density_plot_qc_2 + density_plot_qc_3 + plot_layout(nrow=1) & theme(legend.position="bottom")
scattern_plot_qc_1 <- scattern_plot_qc(so=so_qc, title="Mitochondrial gene percentage", color=pMt_RNA)
scattern_plot_qc_2 <- scattern_plot_qc(so=so_qc, title="Hemoglobin gene percentage", color=pHb_RNA)
scattern_plot_qc_3 <- scattern_plot_qc(so=so_qc, title="Ribsonmal gene percentage", color=pRb_RNA)
scattern_plot_qc_4 <- scattern_plot_qc(so=so_qc, title="Fine labels", color=label_fine_haemosphere) + scale_color_manual(values=color$label_fine_haemosphere)
scattern_plot_qc_5 <- scattern_plot_qc(so=so_qc, title="Doublet class", color=doublet_class) + scale_color_manual(values=c("doublet"="#132B43", "singlet"="#56B1F7"))
scattern_plot_qc_6 <- scattern_plot_qc(so=so_qc, title="QC class", color=qc_class) + scale_color_manual(values=c("fail"="#132B43", "pass"="#56B1F7"))
options(repr.plot.width=22.5, repr.plot.height=10)
scattern_plot_qc_1 + scattern_plot_qc_2 + scattern_plot_qc_3 + scattern_plot_qc_4 + scattern_plot_qc_5 + scattern_plot_qc_6 + plot_layout(ncol=3) & theme(legend.position="bottom")
box_plot_qc_1 <- box_plot_qc(so=so_qc, y=nCount_RNA, fill=tissue, ylab="UMI [count]", ymin=0)
box_plot_qc_2 <- box_plot_qc(so=so_qc, y=nFeature_RNA, fill=tissue, ylab="Feature [count]", ymin=0)
box_plot_qc_3 <- box_plot_qc(so=so_qc, y=pMt_RNA, fill=tissue, ylab="Mt [%]", ymin=0, ymax=100)
box_plot_qc_4 <- box_plot_qc(so=so_qc, y=pHb_RNA, fill=tissue, ylab="Hb [%]", ymin=0, ymax=100)
box_plot_qc_5 <- box_plot_qc(so=so_qc, y=pRb_RNA, fill=tissue, ylab="Rbl [%]", ymin=0, ymax=100)
options(repr.plot.width=22.5, repr.plot.height=20)
box_plot_qc_1[[1]] + box_plot_qc_1[[2]] + box_plot_qc_2[[1]] + box_plot_qc_2[[2]] +
box_plot_qc_3[[1]] + box_plot_qc_3[[2]] + box_plot_qc_4[[1]] + box_plot_qc_4[[2]] +
box_plot_qc_5[[1]] + box_plot_qc_5[[2]] + plot_spacer() + plot_spacer() + plot_layout(guides="collect", ncol=2) & theme(legend.position="none")
so_qc <- subset(so_qc, subset=qc_class=="pass")
so_qc <- NormalizeData(so_qc, assay="RNA")
so_qc <- ScaleData(so_qc, assay="RNA", verbose=FALSE)
so_qc <- RunPCA(so_qc, npcs=50, assay="RNA", features=rownames(so_qc), verbose=FALSE)
so_qc <- FindNeighbors(so_qc, dims=1:30, k.param=20, verbose=FALSE)
so_qc <- FindClusters(so_qc, verbose=FALSE, resolution=1, algorithm=1, group.singletons=FALSE)
so_qc <- RunUMAP(so_qc, dims=1:50, n.neighbors=30, verbose=FALSE)
options(repr.plot.width=22.5, repr.plot.height=15)
dplot_1 <- dplot(so_qc, reduction="umap", group_by="seurat_clusters", label=TRUE) + ggtitle("Louvain cluster")
dplot_2 <- dplot(so_qc, reduction="umap", group_by="tissue") + scale_color_manual(values=color$tissue) + ggtitle("Tissue")
dplot_3 <- dplot(so_qc, reduction="umap", group_by="treatment") + scale_color_manual(values=color$treatment) + ggtitle("Treatment")
dplot_4 <- dplot(so_qc, reduction="umap", group_by="cc_phase_class") + scale_color_manual(values=color$cc_phase_class) + ggtitle("Cell cycle")
dplot_5 <- dplot(so_qc, reduction="umap", group_by="label_main_haemosphere") + scale_color_manual(values=color$label_main_haemosphere) + ggtitle("Haemosphere label (main)")
dplot_6 <- dplot(so_qc, reduction="umap", group_by="label_fine_haemosphere") + scale_color_manual(values=color$label_fine_haemosphere) + ggtitle("Haemosphere label (fine)")
dplot_7 <- dplot(so_qc, reduction="umap", group_by="sample_rep") + ggtitle("Replicate")
fplot_1 <- fplot(so_qc, reduction="umap", features="Trac") + scale_colour_viridis_c(option="inferno") + ggtitle("Trac")
fplot_2 <- fplot(so_qc, reduction="umap", features="Igkc") + scale_colour_viridis_c(option="inferno") + ggtitle("Igkc")
dplot_comb <- dplot_1 + dplot_2 + dplot_3 + dplot_4 + dplot_5 + dplot_6 + dplot_7 + fplot_1 + fplot_2 + plot_layout(ncol=3)
dplot_comb
Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale. Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.
so_qc@meta.data %>% dplyr::group_by(seurat_clusters, treatment, sample_rep) %>%
dplyr::summarise(n=n()) %>% data.frame() %>%
tidyr::spread(seurat_clusters, n) %>%
kableExtra::kable("html") %>% as.character() %>% IRdisplay::display_html()
`summarise()` has grouped output by 'seurat_clusters', 'treatment'. You can override using the `.groups` argument.
| treatment | sample_rep | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13 | 14 | 15 | 16 | 17 | 18 | 19 | 20 | 21 | 22 | 23 | 24 | 25 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| CpG | Rep1 | 41 | 46 | 96 | 922 | 1020 | 213 | 41 | 159 | 162 | 425 | 332 | 7 | 22 | 511 | 130 | 139 | 75 | 180 | 105 | 64 | 68 | 50 | 13 | 4 | 5 | NA |
| CpG | Rep2 | 1455 | 105 | 105 | 105 | 33 | 837 | 823 | 380 | 243 | 11 | 337 | 3 | 604 | 153 | 310 | 251 | 109 | 10 | 58 | 115 | 4 | 10 | 2 | 3 | 1 | 1 |
| NaCl | Rep1 | 3 | 37 | 347 | 259 | 100 | 5 | 18 | 119 | 347 | 422 | 39 | 434 | 2 | 15 | 68 | 71 | 114 | 250 | 124 | 51 | 92 | 30 | 21 | 7 | 9 | 30 |
| NaCl | Rep2 | 266 | 1411 | 777 | 30 | 16 | 79 | 100 | 259 | 153 | 34 | 32 | 279 | 83 | 1 | 122 | 72 | 217 | 7 | 61 | 45 | 16 | 46 | 67 | 44 | 25 | 2 |
marker_22 <- FindMarkers(so_qc, ident.1="22", logfc.threshold=0.5, min.pct=0.01, only.pos=FALSE)
marker_23 <- FindMarkers(so_qc, ident.1="23", logfc.threshold=0.5, min.pct=0.01, only.pos=FALSE)
marker_25 <- FindMarkers(so_qc, ident.1="25", logfc.threshold=0.5, min.pct=0.01, only.pos=FALSE)
options(repr.plot.width=15, repr.plot.height=5)
cluster_marker <- list(marker_22, marker_23, marker_25)
vp_cluster_marker <- lapply(seq_along(cluster_marker), function(i) vp_dea(cluster_marker[[i]], title=c("Cluster 22", "Cluster 23", "Cluster 25")[i], log2_thold=0.25))
wrap_plots(vp_cluster_marker)
so_qc <- subset(so_qc, subset=seurat_clusters!=22)
so_qc <- subset(so_qc, subset=seurat_clusters!=23)
so_qc <- subset(so_qc, subset=seurat_clusters!=25)
# Store data as RDS
saveRDS(so_qc, so_qc_rds)
# Store data as h5ad
adata <- import("anndata", as="adata", convert=FALSE)
pd <- import("pandas", as="pd", convert=FALSE)
np <- import("numpy", as="np", convert=FALSE)
# Transform dgCMatrix to
X <- GetAssayData(so_qc, assay="RNA", slot="counts") %>% as.matrix() %>% t()
X <- np$array(X, dtype=np$int32)
adata <- adata$AnnData(X=X, obs=so_qc@meta.data)
adata$var_names <- rownames(GetAssayData(so_qc, assay="RNA", slot="counts"))
adata$raw <- adata
adata$write_h5ad(so_qc_h5ad)
None
sessionInfo()
R version 4.1.0 (2021-05-18) Platform: x86_64-conda-linux-gnu (64-bit) Running under: Red Hat Enterprise Linux 8.5 (Ootpa) Matrix products: default BLAS/LAPACK: /home/fdeckert/bin/miniconda3/envs/r.4.1.0-FD20200109SPLENO/lib/libopenblasp-r0.3.15.so locale: [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8 [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8 [7] LC_PAPER=en_US.UTF-8 LC_NAME=C [9] LC_ADDRESS=C LC_TELEPHONE=C [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C attached base packages: [1] stats graphics grDevices utils datasets methods base other attached packages: [1] ggrepel_0.9.1 patchwork_1.1.1 RColorBrewer_1.1-3 [4] IRdisplay_1.1 kableExtra_1.3.4 knitr_1.39 [7] biomaRt_2.50.3 reticulate_1.24 forcats_0.5.1 [10] stringr_1.4.0 dplyr_1.0.9 purrr_0.3.4 [13] readr_2.1.2 tidyr_1.2.0 tibble_3.1.7 [16] ggplot2_3.3.6 tidyverse_1.3.1 SeuratWrappers_0.3.0 [19] scDblFinder_1.8.0 sp_1.4-7 SeuratObject_4.1.0 [22] Seurat_4.1.1 loaded via a namespace (and not attached): [1] rappdirs_0.3.3 pbdZMQ_0.3-7 [3] scattermore_0.8 R.methodsS3_1.8.1 [5] bit64_4.0.5 irlba_2.3.5 [7] DelayedArray_0.20.0 R.utils_2.11.0 [9] data.table_1.14.2 rpart_4.1.16 [11] KEGGREST_1.34.0 RCurl_1.98-1.6 [13] generics_0.1.2 BiocGenerics_0.40.0 [15] ScaledMatrix_1.2.0 cowplot_1.1.1 [17] RSQLite_2.2.13 RANN_2.6.1 [19] future_1.25.0 bit_4.0.4 [21] tzdb_0.3.0 spatstat.data_2.2-0 [23] webshot_0.5.3 xml2_1.3.3 [25] lubridate_1.8.0 httpuv_1.6.5 [27] SummarizedExperiment_1.24.0 assertthat_0.2.1 [29] viridis_0.6.2 xfun_0.30 [31] hms_1.1.1 evaluate_0.15 [33] promises_1.2.0.1 fansi_1.0.3 [35] progress_1.2.2 dbplyr_2.1.1 [37] readxl_1.4.0 igraph_1.3.1 [39] DBI_1.1.2 htmlwidgets_1.5.4 [41] spatstat.geom_2.4-0 stats4_4.1.0 [43] ellipsis_0.3.2 RSpectra_0.16-1 [45] backports_1.4.1 deldir_1.0-6 [47] sparseMatrixStats_1.6.0 MatrixGenerics_1.6.0 [49] vctrs_0.4.1 SingleCellExperiment_1.16.0 [51] Biobase_2.54.0 here_1.0.1 [53] Cairo_1.5-15 remotes_2.4.2 [55] ROCR_1.0-11 abind_1.4-5 [57] cachem_1.0.6 withr_2.5.0 [59] progressr_0.10.0 sctransform_0.3.3 [61] prettyunits_1.1.1 scran_1.22.1 [63] goftest_1.2-3 svglite_2.1.0 [65] cluster_2.1.3 lazyeval_0.2.2 [67] crayon_1.5.1 labeling_0.4.2 [69] edgeR_3.36.0 pkgconfig_2.0.3 [71] GenomeInfoDb_1.30.1 nlme_3.1-157 [73] vipor_0.4.5 rlang_1.0.6 [75] globals_0.14.0 lifecycle_1.0.3 [77] miniUI_0.1.1.1 filelock_1.0.2 [79] BiocFileCache_2.2.1 modelr_0.1.8 [81] rsvd_1.0.5 rprojroot_2.0.3 [83] cellranger_1.1.0 polyclip_1.10-0 [85] matrixStats_0.62.0 lmtest_0.9-40 [87] Matrix_1.4-1 IRkernel_1.3 [89] zoo_1.8-10 reprex_2.0.1 [91] base64enc_0.1-3 beeswarm_0.4.0 [93] ggridges_0.5.3 png_0.1-7 [95] viridisLite_0.4.0 bitops_1.0-7 [97] R.oo_1.24.0 KernSmooth_2.23-20 [99] Biostrings_2.62.0 blob_1.2.3 [101] DelayedMatrixStats_1.16.0 parallelly_1.31.1 [103] spatstat.random_2.2-0 S4Vectors_0.32.4 [105] beachmat_2.10.0 scales_1.2.0 [107] memoise_2.0.1 magrittr_2.0.3 [109] plyr_1.8.7 ica_1.0-2 [111] zlibbioc_1.40.0 compiler_4.1.0 [113] dqrng_0.3.0 fitdistrplus_1.1-8 [115] cli_3.4.1 XVector_0.34.0 [117] listenv_0.8.0 pbapply_1.5-0 [119] MASS_7.3-57 mgcv_1.8-40 [121] tidyselect_1.1.2 stringi_1.7.6 [123] highr_0.9 BiocSingular_1.10.0 [125] locfit_1.5-9.5 grid_4.1.0 [127] tools_4.1.0 future.apply_1.9.0 [129] parallel_4.1.0 rstudioapi_0.13 [131] uuid_1.1-0 bluster_1.4.0 [133] metapod_1.2.0 gridExtra_2.3 [135] farver_2.1.0 Rtsne_0.16 [137] digest_0.6.29 BiocManager_1.30.17 [139] rgeos_0.5-9 shiny_1.7.1 [141] Rcpp_1.0.8.3 GenomicRanges_1.46.1 [143] broom_0.8.0 scuttle_1.4.0 [145] later_1.3.0 RcppAnnoy_0.0.19 [147] httr_1.4.3 AnnotationDbi_1.56.2 [149] colorspace_2.0-3 rvest_1.0.2 [151] XML_3.99-0.9 fs_1.5.2 [153] tensor_1.5 IRanges_2.28.0 [155] splines_4.1.0 uwot_0.1.11 [157] statmod_1.4.36 spatstat.utils_2.3-0 [159] scater_1.22.0 xgboost_1.6.0.1 [161] plotly_4.10.0 systemfonts_1.0.4 [163] xtable_1.8-4 jsonlite_1.8.0 [165] R6_2.5.1 pillar_1.8.1 [167] htmltools_0.5.2 mime_0.12 [169] glue_1.6.2 fastmap_1.1.0 [171] BiocParallel_1.28.3 BiocNeighbors_1.12.0 [173] codetools_0.2-18 utf8_1.2.2 [175] lattice_0.20-45 spatstat.sparse_2.1-1 [177] curl_4.3.2 ggbeeswarm_0.6.0 [179] leiden_0.3.10 survival_3.3-1 [181] limma_3.50.3 rmarkdown_2.14 [183] repr_1.1.4 munsell_0.5.0 [185] GenomeInfoDbData_1.2.7 haven_2.5.0 [187] reshape2_1.4.4 gtable_0.3.0 [189] spatstat.core_2.4-2